      SUBROUTINE SHADTA(ELM,SUNPOS,KUMB,KPEN,UMBIN,UMBOUT,PENIN,
     1     PENOUT,IERPNT,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C
C  THIS ROUTINE FINDS THE TRUE ANOMALIES FOR THE ENTRANCE/EXIT TO
C  UMBRA/PENUMBRA SHADOW BY A SPACECRAFT IN ORBIT ABOUT A CENTRAL
C  BODY.
C
C  THE DATA STATEMENTS BELOW DEFINE RC, THE RADIUS OF THE CENTRAL BODY,
C  AS THAT OF THE EARTH. OTHER CENTRAL BODIES CAN BE USED BY CHANGING
C  THE VALUE OF RC.
C
C
C  VAR    DIM  TYPE  I/O  DESCRIPTION
C  ---    ---  ----  ---  -----------
C
C  ELM     6   R*8    I   THE KEPLERIAN ELEMENTS DESCRIBING THE 
C                         SPACECRAFT ORBIT ABOUT THE CENTRAL BODY. 
C                         ORDER IS SMA, ECC, INCL, NODE, ARGP, AND 
C                         ELM(6). ELM(6) IS NOT USED, SO IT DOES NOT 
C                         MATTER WHAT IS THERE. IN INERTIAL COORDINATES.
C                         UNITS ARE KM AND RADIANS.
C
C  SUNPOS  3   R*8    I   POSITION OF THE SUN IN THE SAME COORDINATE
C                         SYSTEM AS ELM. THE POSITION OF THE SUN
C                         SHOULD CORRESPOND TO THE TIME ASSOCIATED WITH
C                         ELM. IN KILOMETERS.
C
C  KUMB    1   I*4   I/O  INPUT/OUTPUT FLAG FOR UMBRA(DEEPEST SHADOW)
C                         CALCULATIONS.
C
C                     I   = NEGATIVE MEANS DO NOT DO UMBRA CALCULATIONS.
C                           IF NEGATIVE, THEN THE VARIABLE IS NOT
C                           RESET. IN THIS CASE THE CALLING PROGRAM
C                           CAN USE AN INTEGER CONSTANT IN THE CALL.
C
C                           IF ZERO OR POSITIVE, LOOK FOR UMBRA ENTRY
C                           AND EXIT TRUE ANOMALIES.
C                           IF NOT NEGATIVE, THEN KUMB RETURNS WITH
C                           ITS VALUE SET TO 0, 1, OR 2.
C
C                     O   = 0 MEANS S/C ESCAPES CENTRAL BODY, IMPACTS
C                             IT, OR AN ALGORITHM ERROR HAS OCCURRED.
C                             NO USEFUL TRUE ANOMALIES ARE RETURNED.
C
C                     O   = 1 MEANS NO SHADOW PASSAGE OCCURS.
C
C                     O   = 2 MEANS SHADOW PASSAGES OCCUR AND THE
C                             UMBRA ENTRY AND EXIT TRUE ANOMALIES ARE
C                             IN UMBIN AND UMBOUT.
C
C  KPEN    1   I*4   I/O  SAME AS FOR KUMB, EXCEPT USED FOR PENUMBRA
C                         (PARTIAL SHADOW).
C
C  UMBIN   1   R*8    O   IF KUMB=2, UMBIN IS THE TRUE ANOMALY OF THE
C                         S/C WHEN IT ENTERS UBMRA(DEEPEST SHADOW).
C                         OTHERWISE, NO USEFUL INFORMATION.
C                         RADIANS.
C
C  UMBOUT  1   R*8    O   IF KUMB=2, UMBOUT IS THE TRUE ANOMALY OF THE
C                         S/C WHEN IT LEAVES UMBRA.
C                         OTHERWISE, NO USEFUL INFORMATION.
C                         RADIANS.
C
C  PENIN   1   R*8    O   IF KPEN=2, PENIN IS THE TRUE ANOMALY OF THE
C                         S/C WHEN IT ENTERS PENUMRA(PARTIAL SHADOW).
C                         OTHERWISE, NO USEFUL INFORMATION.
C                         RADIANS.
C
C  PENOUT  1   R*8    O   IF KPEN=2, PENOUT IS THE TRUE ANOMALY OF THE
C                         S/C WHEN IT LEAVES PENUMBRA.
C                         OTHERWISE, NO USEFUL INFORMATION.
C                         RADIANS.
C
C  IERPNT  1   I*4    I   FORTRAN UNIT NUMBER FOR ERROR MESSAGE PRINTOUT
C                         IN CASE AN ALGORITHM ERROR IS SENSED. IF ZERO,
C                         NO MESSAGE IS GIVEN.
C
C  IERR    1   I*4    O   ERROR RETURN FLAG. 
C                         =0, NO ERROR. =1, ERROR FOUND.
C
C***********************************************************************
C
C  WRITTEN BY RON COOK OF CSC. GIVEN TO CJP 7/77.
C  TRANSFERRED TO VAX WITH MINOR MODS. C PETRUZZO. 9/81.
C
C***********************************************************************
C
C
C THE TECHNIQUE IS PROBABLY SIMILAR TO THAT FOUND IN THE BOOK
C   METHODS OF ORBIT DETERMINATION -- P.R.ESCOBAL,1965, PAGES 155-162.
C
      DIMENSION ELM(6),SUNPOS(3),VEC(3),COEF(5),ROOTS(4),DUMROT(3,3)
      DATA ETOL/0.00150D0/
      REAL*8 RC/6378.14D0/
      REAL*8 RSUN/695500.D0/
      REAL*8 PI/ 3.141592653589793D0 /
      REAL*8 DEGRAD/ 57.29577951308232D0 /
      REAL*8 VCHECK(3),PERIGE(3)
C
      P=ELM(1)*(1.D0-ELM(2)**2)
      E = ELM(2)
      KKUMB=0
      IF(KUMB.GE.0) KKUMB=1
      KKPEN=0
      IF(KPEN.GE.0) KKPEN=1
      TA2IN = 0.D0
      TA2OUT = 0.D0
      TA1IN = 0.D0
      TA1OUT = 0.D0
      IERR=0
      INFO=1
      KOUNT = 0
      ROMIN = P/(1.D0+E)
      IF (ROMIN .GT. RC .AND. E .LT. 1.D0) GO TO 10
C     VEHICLE IMPACTS OR ESCAPES CENTRAL PLANET
      JUMB = 1
      JPEN = 1
      INFO=0
      GO TO 999
C
   10 CONTINUE
C  INITIALIZE: EXPRESS SUN IN COORDS RELATIVE TO ORBIT PLANE.
C
      AINCL=ELM(3)
      ANODE=ELM(4)
      AARGP=ELM(5)
      CW = DCOS(AARGP)
      SW = DSIN(AARGP)
      CI = DCOS(AINCL)
      SI = DSIN(AINCL)
      CO = DCOS(ANODE)
      SO = DSIN(ANODE)
C
C     FIRST, SEE IF INERTIAL Z AXIS, SUN VECTOR, AND PERIGEE VECTOR ARE
C     COPLANAR. IF SO, ALTER ARGUMENT OF PERIGEE SLIGHTLY TO MAKE THEM
C     NON-COPLANAR. EQUATIONS FAIL FOR COPLANAR CASE, HENCE MODS NEEDED.
      PERIGE(1)=CW*CO - SW*SO*CI
      PERIGE(2)=CW*SO + SW*CO*CI
C     PERIGE(3)=SW*SI
C     SEE IF THE PROJECTIONS OF THE SUN AND PERIGEE VECTORS ARE
C     COLLINEAR. IF SO, THEN VECTORS ARE COPLANAR.
      TEMP=DSQRT((PERIGE(1)**2+PERIGE(2)**2) *
     1              (SUNPOS(1)**2+SUNPOS(2)**2))
      IF(TEMP.NE.0) TEMP=(PERIGE(1)*SUNPOS(1)+PERIGE(2)*SUNPOS(2))/TEMP
      IF(DABS(TEMP).LT.(1.D0-1.D-6)) GO TO 9
      DARGP=0.0001D0/DEGRAD
      CW=DCOS(AARGP+DARGP)
      SW=DSIN(AARGP+DARGP)
    9 CONTINUE
      DUMROT(1,1) = CW*CO - SO*SW*CI
      DUMROT(1,2) =CW*SO + CO*SW*CI
      DUMROT(1,3) = SW*SI
      DUMROT(2,1) =-SW*CO-SO*CW*CI
      DUMROT(2,2) =-SW*SO + CO*CW*CI
      DUMROT(2,3) = CW*SI
      DUMROT(3,1) = SI*SO
      DUMROT(3,2) =-SI*CO
      DUMROT(3,3) = CI
      DO 8822 I=1,3
      VEC(I)=DUMROT(I,1)*SUNPOS(1)+DUMROT(I,2)*SUNPOS(2)
     1         +DUMROT(I,3)*SUNPOS(3)
 8822 VEC(I)=-VEC(I)
      CALL UNIVEC(VEC)
      CD = DSQRT(VEC(1)**2+VEC(2)**2)
      TEMP=1.D0-CD*CD
      IF(TEMP.LT.0.D0) TEMP=0.D0
      SD = DSQRT(TEMP)
      CG = VEC(1)/CD
      SG = VEC(2)/CD
      GAMA = DATAN2(VEC(2),VEC(1))
C
C     INITIALIZE FOR UMBRA
      DIS=DSQRT( SUNPOS(1)**2 + SUNPOS(2)**2 + SUNPOS(3)**2 )
      SA =(RSUN - RC)/DIS
      CA = DSQRT(1.D0-SA*SA)
      ALFA=DASIN(SA)
   11 IF (KKUMB .EQ. 0) GO TO 99
    1 IF (KOUNT . NE. 0 .AND. KKPEN.EQ.0) GO TO 99
      IF(KKUMB .NE. 0 .AND.KOUNT.EQ.1) GO TO 12
      IF(KOUNT.EQ.0) GO TO 2
C     CHANGE FOR PENUMBRA
   12 SA = (RSUN + RC)/DIS
      CA = DSQRT(1.D0-SA*SA)
      ALFA = DASIN(SA)
    2 CONTINUE
C     CHECK FOR NO POSSIBLE SHADOW
C     DETERMINE MAX DISTANCE TO SHADOW
      RSMAX = (SA*CD + SD*CA)
      IF(KOUNT .EQ. 1) RSMAX = (SD*CA - SA*CD)
      IF(RSMAX .GT. 1.D-15) GO TO 21
C     MAXIMUM PENUMBRA DISTANCE IS INFINITE
      RSMAX = 1.D20
      GO TO 22
   21 RSMAX = RC / RSMAX
      ROMIN = P/(1.D0+E)
      IF(RSMAX .GT. ROMIN) GO TO 22
C     NO SHADOW POSSIBLE
      IF (KOUNT .EQ. 0) JUMB = 2
      IF(KOUNT .EQ.1)JPEN = 2
      GO TO 99
   22 IF (E .LT. ETOL) GO TO 993
C
C     SET UP COEFFICIENTS FOR QUARTIC IN COS(SI)
C
      IF(KOUNT .EQ. 0) Q=1.0
      IF(KOUNT .EQ. 1) Q=-1.0
      A = (P*Q*SA*CD-RC*E*CG)**2+(RC*E*SG)**2+(P*CA*CD)**2
      B = 2.D0*RC*(P*Q*SA*CD-RC*E*CG)
      C = RC*RC*(1.D0-(E*SG)**2)-(P*CA)**2
      D = 2.D0*P*RC*E*SG*CA
      COEF(1) = C*C-D*D
      COEF(2) = -2.D0*B*C
      COEF(3) = 2.D0*A*C+B*B+D*D*(1.D0+CD*CD)
      COEF(4) = -2.D0*A*B
      COEF(5) = A*A - (D*CD)**2
      CALL ROOT34(COEF,4,ROOTS,NROOTS)
      NUM = NROOTS + 1
      GO TO (4,4,5,4,5),NUM
    4 CONTINUE
C     NO REAL ROOTS(ONE OR THREE NOT POSSIBLE)
      IF(KOUNT .EQ. 0) JUMB = 5
      IF(KOUNT .EQ. 1) JPEN = 5
      GO TO 99
    5 CONTINUE
C     TWO OR FOUR REAL ROOTS
C     CHECK QUADRANTS OF SOLUTIONS
      DO 51 J=1,NROOTS
      RADICL=1.D0-ROOTS(J)**2
      IF(RADICL.LT.0.D0) GO TO 7630
      SSI=DSQRT(RADICL)
      CBETA = CD*ROOTS(J)
      SBETA = DSQRT(1.D0-CBETA**2)
      RPLS = P/(1.D0+E*(ROOTS(J)*CG-SSI*SG))
      RMNS = P/(1.D0+E*(ROOTS(J)*CG+SSI*SG))
      RSHAD = RC/(SA*CBETA*Q+CA*SBETA)
      DELRP = DABS(RPLS-RSHAD)
      DELRM = DABS(RMNS - RSHAD)
      IF(DELRM .LT. DELRP) SSI=-SSI
C     REPLACE COS(SI) WITH SI
   51 ROOTS(J) = DATAN2(SSI,ROOTS(J))
 8888 CONTINUE
      IF (KOUNT .EQ. 1) GO TO 6
C
C     SELECT ANGLES FOR UMBRAL BOUNDARIES
C
      DO 52 J=1,NROOTS
      DO 52 K=J,NROOTS
      IF(DABS(ROOTS(J)).LT.DABS(ROOTS(K))) GO TO 52
      TEMP = ROOTS(J)
      ROOTS(J) = ROOTS(K)
      ROOTS(K) = TEMP
   52 CONTINUE
C     CHECK THAT SOLUTIONS ARE ON NITE SIDE OF UMBRAL TERMINATOR
      TEST = 0.5D0*PI-ALFA
      IF(DABS(ROOTS(1)).LT.TEST) GO TO 53
C     NO UMBRAL INTERSECTIONS
      JUMB = 3
      GO TO 99
   53 TA2IN = ROOTS(1)
      TA2OUT = ROOTS(2)
      IF(TA2OUT .GT. TA2IN) GO TO 54
      TA2IN = ROOTS(2)
      TA2OUT = ROOTS(1)
   54 JUMB = 4
      GO TO 99
C     SELECT ANGLES FOR PENUMBRAL BOUNDARIES
    6 DO 61 J=1,NROOTS
      DO 61 K=J,NROOTS
      IF (DABS(ROOTS(J)).LT.DABS(ROOTS(K))) GO TO 61
      TEMP = ROOTS(J)
      ROOTS(J) = ROOTS(K)
      ROOTS(K) = TEMP
   61 CONTINUE
C     CHECK TO MAKE SURE THAT SOLUTIONS ARE ON NITE SIDE
C     OF PENUMBRAL TERMINATOR
      TEST = 0.5D0*PI+ALFA
      IF(DABS(ROOTS(1)).LT.TEST) GO TO 62
C     NO PENUMBRAL INTERSECTIONS
      JPEN = 3
      GO TO 99
   62 TA1IN = ROOTS(1)
      TA1OUT = ROOTS(2)
      IF(TA1IN.LT.TA1OUT) GO TO 63
      TA1IN = ROOTS(2)
      TA1OUT = ROOTS(1)
   63 JPEN = 4
   99 KOUNT = KOUNT + 1
      IF(JUMB .NE. 4 .OR. KOUNT .NE. 1) GO TO 991
      TA2IN = TA2IN + GAMA
      TA2OUT = TA2OUT + GAMA
  991 IF (JPEN.NE.4.OR.KOUNT.NE.2) GO TO 992
      TA1IN = TA1IN + GAMA
      TA1OUT = TA1OUT + GAMA
  992 IF(KOUNT .EQ. 2) GO TO 999
      IF(KOUNT .EQ. 1 .AND. KKPEN .NE. 0) GO TO 1
      GO TO 999
C     APPROXIMATE SOLUTION FOR NEAR-CIRCULAR ORBITS
  993 CBETA = RC*SA/ROMIN
      SBETA = CA*DSQRT(1.D0-(RC/ROMIN)**2)
      IF (KOUNT .NE. 0) GO TO 994
      CBETA = CBETA +SBETA
C     UMBRA
      SI = DACOS(CBETA/CD)
      TA2IN = -SI
      TA2OUT = SI
      JUMB = 4
      GO TO 99
  994 CBETA = SBETA - CBETA
C     PENUMBRA
      SI = DACOS(CBETA/CD)
      TA1IN = -SI
      TA1OUT = SI
      JPEN = 4
      GO TO 99
C
 7630 CONTINUE
      IERR=1
      INFO=0
      IF(IERPNT.NE.0) WRITE(IERPNT,7632)
 7632 FORMAT(' SHADTA. SHADOW ALGORITHM ERROR.')
C
C
  999 CONTINUE
C
C     NOW SET UP UMBRA/PENUMBRA RETURN FLAGS.
C
C     JUMB/JPEN FLAGS HAVE BEEN DEFINED AS FOLLOWS ----
C
C     JUMB  UBMRA(DEEP SHADOW) FLAG. 
C     JPEN  PENUMBRA(PARTIAL SHADOW) FLAG. 
C
C     JUMB/JPEN ASSIGNMENTS HAVE BEEN --
C     =1  S/C IMPACTS OR ESCAPES PLANET. NO SHADOW
C         INFO RETURNED.
C     =2  NO SHADOW. S/C     DOES NOT ENTER SHADOW.
C     =3  NO REAL SOLUTIONS ON NIGHT SIDE OF
C         TERMINATOR. EFFECTIVELY, NO SHADOW
C         PASSAGE OCCURS. 
C     =4  ENTRY AND EXIT TRUE ANOMALIES FOUND.
C         SHADOW OCCURS.
C     =5  NO REAL SOLUTIONS. EFFECTIVELY, NO SHADOW
C         PASSAGE OCCURS.
C
C     DO PENUMBRA RETURN FLAG.
      IF(KPEN.LT.0) GO TO 800
      PENIN=999.D0/DEGRAD
      PENOUT=PENIN
      IF(INFO.NE.0) GO TO 710
C     ERROR, ESCAPE, OR IMPACT. NO INFO.
      KPEN=0
      GO TO 800
  710 CONTINUE
      KPEN=1
      IF(JPEN.EQ.4) KPEN=2
      IF(KPEN.NE.2) GO TO 800
      PENIN=TA1IN
      PENOUT=TA1OUT
C 
  800 CONTINUE
C     DO UMBRA RETURN FLAG.
      IF(KUMB.LT.0) GO TO 900
      UMBIN=999.D0/DEGRAD
      UMBOUT=UMBIN
      IF(INFO.NE.0) GO TO 810
C     ERROR, ESCAPE, OR IMPACT. NO INFO.
      KUMB=0
      GO TO 900
  810 CONTINUE
      KUMB=1
      IF(JUMB.EQ.4) KUMB=2
      IF(KUMB.NE.2) GO TO 900
      UMBIN=TA2IN
      UMBOUT=TA2OUT
C
  900 CONTINUE
      RETURN
      END
